{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dcb68ef2-d318-4f51-bfdf-bda799bf9efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import cantera as ct\n",
    "import matplotlib\n",
    "import scipy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "42b576f0-ab50-4e69-8673-a7cbcf656b92",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameter values\n",
    "p0 = 101325  # pressure\n",
    "T0 = 296  # inlet temperature\n",
    "phi=0.7 #equivalence ratio\n",
    "alpha=0.566 #<1 volume fraction\n",
    "mdot_reactants = 80  # kg/m^2/s\n",
    "mdot_products = mdot_reactants  # kg/m^2/s\n",
    "\n",
    "# Mechanism and fuel composition\n",
    "rxnmech = 'gri30.yaml'  # reaction mechanism file\n",
    "xh2 = alpha / (1 - alpha) * 1  # xh2 = alpha/(1-alpha)*xc\n",
    "fuel = {'CH4': 1, 'H2': xh2}\n",
    "oxidizer = {'O2': 1, 'N2': 3.76}\n",
    "\n",
    "# Calculation params\n",
    "width = 0.1  # m domain width\n",
    "loglevel = 0  # amount of diagnostic output (0 to 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1a4df806-0aea-4c9e-93e4-fc701806c3f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gas object 1 for the counterflow flame\n",
    "gas1=ct.Solution(rxnmech)\n",
    "gas1.TP = T0, p0\n",
    "#fix the composition of the fuel\n",
    "gas1.set_equivalence_ratio(phi, fuel, oxidizer)  # hold temperature and pressure constant \n",
    "\n",
    "# Create the counterflow premixed flame simulation object\n",
    "fl1 = ct.CounterflowPremixedFlame(gas=gas1, width=width)\n",
    "fl1.transport_model = 'multicomponent'\n",
    "fl1.energy_enabled = True #energy equation\n",
    "fl1.set_refine_criteria(ratio=3, slope=0.1, curve=0.2, prune=0.02)\n",
    "\n",
    "# set the boundary flow rates\n",
    "fl1.reactants.mdot = mdot_reactants\n",
    "fl1.products.mdot = mdot_products\n",
    "\n",
    "fl1.set_initial_guess()  # assume adiabatic equilibrium products\n",
    "fl1.solve(loglevel, auto=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1b933a4c-484e-4e13-bb0a-cace74025608",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2000.0\n",
      "2000.0\n"
     ]
    }
   ],
   "source": [
    "# Gas object 2 for the counterflow flame with non equilibrium products\n",
    "gas2=ct.Solution(rxnmech)\n",
    "gas2.TP = T0, p0\n",
    "gas2.set_equivalence_ratio(phi, fuel, oxidizer)  # hold temperature and pressure constant\n",
    "\n",
    "# Create the flame simulation object\n",
    "fl2 = ct.CounterflowPremixedFlame(gas=gas2, width=width)\n",
    "# Set grid refinement parameters\n",
    "fl2.set_refine_criteria(ratio=3, slope=0.1, curve=0.2, prune=0.02)\n",
    "# set the boundary flow rates\n",
    "fl2.reactants.mdot = mdot_reactants\n",
    "fl2.products.mdot = mdot_products\n",
    "fl2.products.T =2000 #product temperature, temperature of my project\n",
    "fl2.products.X =fl1.products.X #composition of the products\n",
    "fl2.set_initial_guess(equilibrate=False)\n",
    "print(fl2.products.T)\n",
    "fl2.solve(loglevel, auto=True)\n",
    "print(fl2.products.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4bcdc4d-ad6a-493a-b750-32f81a61c4c7",
   "metadata": {},
   "source": [
    "## Simple example using counterflow premixed flame in Cantera"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f09d899d-0634-43fa-8ed8-fdcde5162a31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAGdCAYAAADT1TPdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/lUlEQVR4nO3dd3xUVf7/8ffMJJn0EEhIgQBJ6B2D0kVXF1TEthYsrPgFbLCi4vqDtYAF0EVXV1dFUVFxV2xrQRErKlJEkEivCSQQQkiATAppM/f3RyBLTIAEZnJnktfz8ZiHZM6dO597DDNv7j3nXIthGIYAAAC8lNXsAgAAAE6GsAIAALwaYQUAAHg1wgoAAPBqhBUAAODVCCsAAMCrEVYAAIBXI6wAAACv5md2AWfK5XIpKytLYWFhslgsZpcDAADqwDAMFRQUKD4+Xlbryc+d+HxYycrKUkJCgtllAACA05CZmanWrVufdBufDythYWGSKg82PDzc5GoAAEBdOBwOJSQkVH2Pn4zPh5Vjl37Cw8MJKwAA+Ji6DOFggC0AAPBqhBUAAODVCCsAAMCrEVYAAIBXI6wAAACvZnpYmT59uiwWS7VHbGys2WUBAAAv4RVTl7t166Zvvvmm6mebzWZiNQAAwJt4RVjx8/PjbAoAAKiV6ZeBJGn79u2Kj49XYmKiRo0apbS0tBNuW1paKofDUe0BAAAaL9PDSr9+/fTWW2/pyy+/1Ny5c5Wdna2BAwcqLy+v1u1nzZqliIiIqgf3BQIAoHGzGIZhmF3E8YqKipScnKz7779f9957b4320tJSlZaWVv187N4C+fn5LLcPAICPcDgcioiIqNP3t1eMWTleSEiIevTooe3bt9fabrfbZbfbG7gqAABgFq8LK6Wlpdq8ebOGDBliah2/ZhzSJ2v3KsTuV/kIsCnE7qdQu5+C7X4KtduOPn+03W6T3Y9ZTAAAuJvpYeW+++7TyJEj1aZNG+Xk5Ojxxx+Xw+HQzTffbGpdm/c59OaK3fV6jb/NclyAqQwzYYH+SogMUmJUiJKjQ5UYFaLWkUHys5k+XAgAAJ9geljZs2ePrr/+euXm5io6Olr9+/fXypUr1bZtW1Pr6hYfob/8ob0KSytUVFqhojJn5X9LK1RYWvnn4rIKFZZWqKTcJUkqdxo6XFyuw8XlJ923v82iNs2DlRgVqqToECVFhSgxKkRJ0aGKCg2o0+2yAQBoKrxugG191WeAjqdUOF1VYaYywDiPhpoK5R8pV0ZesdJyC5V2oEi78oqqwk1twux+SqwKMKFKjA5RcnSIOsWEcTYGANBo+PQAW1/kZ7MqIsiqiCD/U27rchna5yhR+oGiqgCTlluk9NxC7Tl0RAWlFVq3J1/r9uRXe12o3U9nt4tU/6QW6p/UQt3iwwkvAIAmgTMrXqSk3KmMg8VHA0zh0UBTpG37C1RQUlFtW8ILAMCX1ef7m7DiA5wuQ5v3ObQyLU8r0w5qVXqeHIQXAIAPI6w0coQXAICvI6w0MXUJL2F2Pw3rFqvLe8drYHILggsAwFSElSbuVOElKjRAI3rE6bLe8TqrTSRTpQEADY6wgmqcLkNrdh/Sp7/t1aL12TpYVFbV1joySJf1itdlvePVOZb+AwA0DMIKTqjc6dJPO3K1MDVLX27MVlGZs6qtU0yYLusdr8t6xSuhebCJVQIAGjvCCurkSJlT323J0Sepe/X91gMqc/5vsbqz2jTTZb3iNaJnvKLDuHEkAMC9CCuot/wj5fpyQ7Y++W2vVuzMk+vob4XVIg1qH6XLesVrePdYhQeeeuE7AABOhbCCM5LjKNFn6/bp09+ylJp5uOp5u59Vf0pprbGDE5UcHWpegQAAn0dYgdvszivSwt+y9ElqlrbnFEqSLBbpgs4xGj8kUeckNmc2EQCg3ggrcDvDMLQq/aDmLk3Xt1v269hvTc/WERo3JEmXdI9l7RYAQJ0RVuBROw8U6rWf0vXhmj0qragclNuqWZBuGdRO152doDDGtQAAToGwggaRV1iqt1dm6K0Vu5R3dO2WMLufru/XRmMGtlN8syCTKwQAeCvCChpUSblTH6/dq7lL07TzQJEkyc9q0YiecRo/JEndW0WYXCEAwNsQVmAKl8vQ99tyNPfHdK1Iy6t6vn9Sc40fkqTzO7WU1cpgXAAAYQVeYMPefM1dmqbP1u2T8+iiLcnRIRo/JEl/SmktfwbjAkCTRliB18g6fERvLN+ld37OUEFp5c0U2zQP1r1/7KiRveJl40wLADRJhBV4nYKSci1YlamXf9yp3MLKwbidYsJ03/BOurBLS9ZqAYAmhrACr1VcVqF5y3Zpzg87VVBSeaald0Iz3T+8kwa2jzK5OgBAQyGswOvlF5fr5R93at6yXTpSXnnn58Hto3Tf8E7qndDM3OIAAB5HWIHPyCko0Qvf7dB/VmWo3Fn5qziiR5ymXNxZCc2DTa4OAOAphBX4nMyDxXr2m+36aO0euQwpwM+qW4ck6Y7zkhVi9zO7PACAmxFW4LM273Po0YWbqtZpiQm36/7hnXVln1as0QIAjQhhBT7NMAx9uXG/Zi7arIyDxZKkXgnN9PClXZXSNtLt75edX6JHFm5U1uEjGtqppcYPSeT+RgDgYYQVNAol5U69vixdL3y3Q0VllYNwr+gdr/93cWfFRbjvvkO3zV+tLzfur/o5PiJQL92Uol4M9AUAj6nP9zfLiMJrBfrbdOd57bXkvvN0TUprWSzSx6lZ+sNTP2juj2mqcLrO+D3KKlz6bkuOJOmeCzuqTfNgZeWX6Pq5K7V0+4Ez3j8A4MwRVuD1WoYHavY1vfTphMHq2zZSR8qdmrFosy771zL9lnn4jPaddfiIyp2GAv2tuuuC9lo0aYgGt49ScZlTY99crVXpB91zEACA00ZYgc/o0TpC7902QE/+qYcigvy1aZ9DV7y4TNM/3aiCkvLT2uexMTFtmgfLYrEo1O6n18b01YVdYlRW4dK4N3/R1uwCdx4GAKCeCCvwKVarRded3UbfTh6qK/u0kmFIbyzfpQv/8YMWb9in+g7BOj6sHGP3s+lfN/TR2e0i5Sip0Pi3Viu/+PTCEADgzBFW4JOiQu165rreentsP7VrEaz9jlLd/vavGv/Wau09fKTO+8k8GlZ+vwBdoL9Nc//cV60jg5RxsFj3vpcql8unx6IDgM8irMCnDe4QpcV3n6u//KG9/G0WfbM5R3/8xw+atyy9TuGitjMrxzQLDtCcm1IU4GfVt1ty9NIPO91ePwDg1Agr8HmB/jZNHtZJX0waonPaNVdxmVOPLNyk6+euVEZe8Ulfe7KwIkndW0Xo8cu7S5Ke+Xqb1u057NbaAQCnRlhBo9G+ZZgW3Npfj13RXcEBNv2cflAX/fNHzV+5u9azLIZhVIWZE4UVSbqmb2uN6BGnCpehe95NVcnRGy8CABoGYQWNitVq0ej+bbV40rnql1h5luWhjzdo9Os/a8+h6mdZDhWXq6C0QlLNMSvHs1gsevyK7moZZtfOA0V6cvEWjx4DAKA6wgoapTYtgvXO+P6aPrKrAv2tWrYjT8Of+VHvrMqomjG0O69IkhQbHqhAf9tJ9xcZEqC/X91TUuXso9QzXN8FAFB3hBU0WlarRWMGJWrxpHPVt22kisqcmvrf9frz66u051CxtucUSpLatjjxWZXjndepZdV06b/9d71bVtAFAJwaYQWNXruoEL172wA9OKKL7H5WLd2eq+HP/KgHP94gSfW6B9ADI7pULUj3xvJdnikYAFANYQVNgs1q0bghSfpi0hCd3a7yLEtZReWZkaEdo+u8n6hQu6Ze3FmS9I+vt9VrTRcAwOnxirDy4osvKjExUYGBgUpJSdHSpUvNLgmNVFJ0qN69dYAevbybusWHa8L5yRqY3KJe+7i2b4L6to1UcZlTj3+2yUOVAgCOMT2svPvuu7r77rv1wAMPaO3atRoyZIguvvhiZWRkmF0aGimr1aI/D2inz+8aor8O7yyLxVLv18+4soesFumLDdlamZbnoUoBAJIXhJV//OMfGjt2rMaNG6cuXbro2WefVUJCgl566SWzSwNOqFNsmG7o10aS9OjCTXKyFD8AeIyfmW9eVlamNWvWaMqUKdWeHzZsmJYvX17ra0pLS1VaWlr1s8PhkCSlpqYqNDS06vnIyEglJiaqpKREmzbVPFV/1llnSZK2bt2qoqKiam3t2rVT8+bNdeDAAWVmZlZrCwsLU4cOHeR0OvXbb7/V2G+PHj3k7++vnTt3Kj8/v1pbq1atFBMTo0OHDik9Pb1aW1BQkLp06SJJWrt2bY0b8nXp0kVBQUHavXu38vKq/0s+JiZGrVq1UkFBgbZv316tzd/fXz169JAkrV+/XuXl1W/I16FDB4WFhWnv3r3av39/tbYWLVqobdu2OnLkiDZv3lytzWKxqE+fPpKkzZs368iR6mM3EhMTFRkZqf3792vv3r3V2iIiIpScnKzy8nKtX79ev9erVy/ZbDZt375dBQXV73ickJCg6OhoHTx4ULt27arWFhISok6dOkmSfv311xr77dq1qwIDA5Wenq5Dhw5Va4uLi1NcXJwcDod27NhRrc1ut6tbt26SpHXr1qmionJtlguiy7XgULo2lMXqgzWZGhRnVU5OTrXXRkVFqU2bNiouLtaWLdXXZ7Farerdu7ckadOmTSopKanWnpSUpGbNmik7O1tZWVnV2po1a6akpCSVlZVpw4YNNY61d+/eslqt2rZtmwoLC6u1tWnTRlFRUcrNza1xBjM0NFQdO3aUy+VSampqjf12795dAQEBSktL0+HDh6u1xcfHKzY2VocPH1ZaWlq1tsDAQHXt2lVS5d9Vl6v6TKrOnTsrODhYGRkZys3NrdbWsmVLtW7dWoWFhdq2bVu1Nj8/P/XsWTmlfOPGjdU+GySpffv2Cg8P1759+7Rv375qbXxGVOIz4n/c/RlxTMeOHRUaGqo9e/bwGXHcZ8Tv3/ekDBPt3bvXkGQsW7as2vMzZswwOnbsWOtrpk2bZkg65ePGG280DMMwtm/fXmv7Mf3796/RNn/+fMMwDONf//pXjbZhw4YZhmEY+fn5te43JyfHMAzDGDlyZI22p59+2jAMw3jvvfdqtPXp06eqpoCAgBrtGzZsMAzDMMaOHVujbcqUKYZhGMaSJUtqtLVq1apqv61atarRvmTJEsMwDGPKlCk12saOHWsYhmFs2LChRltAQEDVfvv06VOj/b333jMMwzCefvrpGm0jR440DMMwcnJyau3D/Px8wzAMY9iwYTXa/vWvfxmGYRjz58+v0da/f/+qmmrb7/bt2w3DMIwbb7yxRtu0adMMwzCMxYsX12hLTk6u2m9UVFSN9tibZhspj31tTPjLpBptd955p2EYhrFmzZoabWFhYVX77dq1a432Tz75xDAMw5g5c2aNtquvvtowDMPIzMys9VhLSkoMwzCMoUOH1mibO3euYRiGMXfu3BptQ4cONQzDMEpKSmrdb2ZmpmEYhnH11VfXaJs5c6ZhGIbxySef1Gjr2rVr1bGGhYXVaF+zZo1hGIZx55131mi75557DMMwjOXLl9doi4qKqtpvcnJyjfbFixcbhlH75wafEXxG/P7hqc+I5cuXG4ZhGPfcc0+NNj4j/vf/82QsR/+nmSIrK0utWrXS8uXLNWDAgKrnZ8yYofnz59dImVLtZ1YSEhL0ww8/cGaFfzU1+L+ayitcuu/rA8oscOn6biG6plt4tdc25X81HY8zK//DZ0SlpvIZIXFm5ZjazqwMHTpU+fn5Cg8Pr7Gv45kaVsrKyhQcHKz3339fV155ZdXzkyZNUmpqqn744YdT7sPhcCgiIqJOBwt4wjeb9mvcW6sVYLPqm3uHqk0dF5kDgKasPt/fpg6wDQgIUEpKir7++utqz3/99dcaOHCgSVUB9XNBl5Ya1L6Fypwu/f1L7hsEAO5m+myge++9V6+++qpef/11bd68Wffcc48yMjJ0++23m10aUCcWi0UPXNJVFov02bp9+jXj0KlfBACoM9PDynXXXadnn31Wjz76qHr37q0ff/xRixYtUtu2bc0uDaizrvHhuvqs1pKkxz/bVGM8AQDg9Jk6ZsUdGLMCb7HfUaLzZn+vI+VOvXDDWRrRM87skgDAa/nMmBWgMYkJD9St5yZJkp5cvEWlFU6TKwKAxoGwArjRrecmKTrMroyDxZq/YrfZ5QBAo0BYAdwoxO6n+4Z1lCQ99+12HSoqM7kiAPB9hBXAza5OSVDn2DA5Sir03HfbT/0CAMBJEVYAN7NZLXpgROVKo/NX7FZ6btEpXgEAOBnCCuABQzpEa2jHaFW4DD35BQvFAcCZIKwAHvLAiC6yWqTFG7O1Kv2g2eUAgM8irAAe0jEmTNed3UaSNOPzTXK5fHpJIwAwDWEF8KB7/thBIQE2/bYnXwvXZZ36BQCAGggrgAe1DAvU7UOTJUl/X7xVJeUsFAcA9UVYATxs3JAkxYYHau/hI5q3bJfZ5QCAzyGsAB4WFGDTX4d3kiS9uGSH8gpLTa4IAHwLYQVoAFf2aaVu8eEqKK3QP79loTgAqA/CCtAArMctFPfvnzO0I6fQ5IoAwHcQVoAGMjA5Shd2aSmny9ATX2w2uxwA8BmEFaABTbm4i2xWi77ZnKPlO3PNLgcAfAJhBWhA7VuG6sZ+lQvFzVy0mYXiAKAOCCtAA5t0QQeF2f20Ya9DH63da3Y5AOD1CCtAA2sRated57eXJM3+cquKyypMrggAvBthBTDBLYPaqXVkkLIdJZrz/U6zywEAr0ZYAUwQ6G/TA5dUTmV++cc07TlUbHJFAOC9CCuASS7qHqv+Sc1VWuHSrC+2mF0OAHgtwgpgEovFoocv7SarRfp83T79nJZndkkA4JUIK4CJusaH6/pzKqcyP7Jwk5xMZQaAGggrgMnu/WNHhQX6adM+h95fnWl2OQDgdQgrgMlahNp194UdJVVOZXaUlJtcEQB4F8IK4AX+PKCtkqNDlFdUpue5KzMAVENYAbyAv82qhy7tKkmat2yX0g5wV2YAOIawAniJ8zq11PmdolXhMjTjc+7KDADHEFYAL/LgpV3lZ7Xo2y05+n5rjtnlAIBXIKwAXiQ5OlRjBraTJD322SaVO13mFgQAXoCwAniZv1zQQS1CArTzQJHmr9htdjkAYDrCCuBlIoL8NXlYJ0nSs99s08GiMpMrAgBzEVYAL3Td2QnqEhcuR0mF/vH1VrPLAQBTEVYAL2SzWjRtZOVU5v/8nKHN+xwmVwQA5iGsAF6qf1ILjegRJ5chPbpwkwyD+wYBaJoIK4AXm3JxZwX4WbUiLU9fbtxvdjkAYArCCuDFEpoH67ZzkyRJMxZtUkm50+SKAKDhEVYAL3fHecmKDQ9U5sEjeuXHNLPLAYAGR1gBvFxwgJ/+NqKLJOmFJTuUebDY5IoAoGGZGlbatWsni8VS7TFlyhQzSwK80sieceqX2FylFS7uGwSgyTH9zMqjjz6qffv2VT0efPBBs0sCvI7FYtEjl3eTzWrR4o3Z+nHbAbNLAoAGY3pYCQsLU2xsbNUjNDTU7JIAr9Q5Nlx/HtBWkjR94UaVVXDfIABNg+lh5cknn1SLFi3Uu3dvzZgxQ2VlLC0OnMjdF3ZUVGiA0g4Uad6ydLPLAYAG4Wfmm0+aNElnnXWWIiMjtWrVKk2dOlXp6el69dVXT/ia0tJSlZaWVv3scLCyJ5qOiCB//b+LOuuvH6zTc99u1+W9Wyk2ItDssgDAo9x+ZmX69Ok1Bs3+/rF69WpJ0j333KOhQ4eqZ8+eGjdunObMmaPXXntNeXl5J9z/rFmzFBERUfVISEhw9yEAXu1PZ7VWnzbNVFTm1KwvGGwLoPGzGG5ewzs3N1e5ubkn3aZdu3YKDKz5r8G9e/eqdevWWrlypfr161fra2s7s5KQkKD8/HyFh4efWfGAj9iwN18j//WTDENacGt/9U9qYXZJAFAvDodDERERdfr+dvtloKioKEVFRZ3Wa9euXStJiouLO+E2drtddrv9tPYPNBbdW0XohnPa6N8/Z+jhTzbo87uGyN9m+hA0APAI0z7dVqxYoWeeeUapqalKT0/Xe++9p9tuu02XXXaZ2rRpY1ZZgM+4b1gnRQb7a9v+QgbbAmjUTAsrdrtd7777rs477zx17dpVDz/8sMaPH6933nnHrJIAnxIZEqCpF1eubPvsN9uVdfiIyRUBgGe4fcxKQ6vPNS+gsXG5DF378gqt3n1Iw7vF6OXRfc0uCQDqpD7f31zkBnyY1WrR41d2l81q0Zcb9+u7LfvNLgkA3I6wAvi4zrHhGjs4UZL08CcbdaTMaXJFAOBehBWgEZh0QQfFRwRqz6EjemHJDrPLAQC3IqwAjUCI3U8Pj+wmSXr5x53akVNockUA4D6EFaCRGN4tRud3ila509BDH2+Qj4+dB4AqhBWgkbBYLHrksu6y+1m1Ii1Pn6RmmV0SALgFYQVoRNq0CNZf/tBekvT455uVf6Tc5IoA4MwRVoBGZvy5SUqKDlFuYame/mqr2eUAwBkjrACNjN3Ppscv7y5Jmr9yt9btOWxuQQBwhggrQCM0sH2ULu8dL8OQHvhog5wuBtsC8F2EFaCRemBEF4UF+mn93nz9++fdZpcDAKeNsAI0Ui3DAvXX4Z0kSbO/3KqcghKTKwKA00NYARqxG/u1Vc/WESooqdDMzzebXQ4AnBbCCtCI2awWPX5Fd1ks0sepWVq+I9fskgCg3ggrQCPXs3Uzje7fVpL04CcbVFrBjQ4B+BbCCtAETB7WSVGhdqUdKNILS3aaXQ4A1AthBWgCIoL89chllTc6fOn7Hdq2v8DkigCg7ggrQBNxSY9YXdilpcqdhqZ8uE4u1l4B4CMIK0ATYbFY9NgV3RVq99OvGYf1NmuvAPARhBWgCYmLCNL9F1WuvfLkF1uUdfiIyRUBwKkRVoAm5qZ+bXVWm2YqKnPq4U82yDC4HATAuxFWgCbGarXoyT/1lL/Nom8252jR+myzSwKAkyKsAE1Qh5gw3Xlee0nStE83Kr+43OSKAODECCtAE3Xn+clq3zJUuYWlmrmIpfgBeC/CCtBE2f1smnVVD0nSu6sztXwnS/ED8E6EFaAJO7tdc93Uv40k6W//Xa+ScpbiB+B9CCtAE3f/RZ0VE27Xrrxi/fPb7WaXAwA1EFaAJi480F+PXd5dkvTKj2nalOUwuSIAqI6wAkDDusXq4u6xcroMTfnvOjlZih+AFyGsAJAkPXJZN4UF+mndnnzNW5ZudjkAUIWwAkCS1DI8UH+7pIsk6emvtinzYLHJFQFAJcIKgCrX9U3QOYnNdaTcqQc+Zil+AN6BsAKgitVq0ayreijAz6oftx3Qf3/da3ZJAEBYAVBdcnSoJl3QQZL0yMKN2u8oMbkiAE0dYQVADbedm6SerSPkKKnQ3/67nstBAExFWAFQg5/Nqqeu6aUAm1XfbsnhchAAUxFWANSqY0yYJl3I5SAA5iOsADghLgcB8AaEFQAnxOUgAN6AsALgpLgcBMBshBUAp8TlIABm8mhYmTFjhgYOHKjg4GA1a9as1m0yMjI0cuRIhYSEKCoqSnfddZfKyso8WRaAeuJyEAAzeTSslJWV6ZprrtEdd9xRa7vT6dSIESNUVFSkn376SQsWLNCHH36oyZMne7IsAKeBy0EAzGIxGuB87htvvKG7775bhw8frvb8F198oUsvvVSZmZmKj4+XJC1YsEBjxoxRTk6OwsPDT7lvh8OhiIgI5efn12l7AKevwunSVS8t17o9+fpD55Z67ea+slgsZpcFwAfV5/vb1DErK1asUPfu3auCiiQNHz5cpaWlWrNmTa2vKS0tlcPhqPYA0DCOvxz0HZeDADQQU8NKdna2YmJiqj0XGRmpgIAAZWdn1/qaWbNmKSIiouqRkJDQEKUCOIrLQQAaWr3DyvTp02WxWE76WL16dZ33V9spZMMwTnhqeerUqcrPz696ZGZm1vcQAJyh42cHTWV2EAAP86vvCyZOnKhRo0addJt27drVaV+xsbH6+eefqz136NAhlZeX1zjjcozdbpfdbq/T/gF4xrHLQZc+91PV5aA/pbQ2uywAjVS9w0pUVJSioqLc8uYDBgzQjBkztG/fPsXFxUmSvvrqK9ntdqWkpLjlPQB4xrHLQbO/3KpHFm7U4A5RigkPNLssAI2QR8esZGRkKDU1VRkZGXI6nUpNTVVqaqoKCwslScOGDVPXrl01evRorV27Vt9++63uu+8+jR8/npk9gA84/nLQXz9YJ5eLy0EA3M+jYeXhhx9Wnz59NG3aNBUWFqpPnz7q06dP1ZgWm82mzz//XIGBgRo0aJCuvfZaXXHFFXrqqac8WRYAN/GzWfX0Nb1k97Pqx20HNG/5LrNLAtAINcg6K57EOiuA+eav2KWHPtmoAJtVH08YpK7x/F0EcHI+s84KgMbhpv5tdWGXlipzujRpwVqVlDvNLglAI0JYAXDGLBaLnvxTT0WH2bU9p1AzPt9sdkkAGhHCCgC3aBFq19PX9JIkzV+5W99u3m9yRQAaC8IKALc5t2O0xg5OlCT99YN1yilgdVsAZ46wAsCt7r+ok7rEhetgUZnue5/pzADOHGEFgFvZ/Wx6blRvpjMDcBvCCgC36xATpgcv7SpJevKLLdqUxd3RAZw+wgoAj7ipXxtd2CWG6cwAzhhhBYBHVE5n7sF0ZgBnjLACwGN+P535m01MZwZQf4QVAB51/HTm+z9cpxwH05kB1A9hBYDHHT+defL7vzGdGUC9EFYAeNzx05mXbs/V68vSzS4JgA8hrABoEMdPZ/774q1MZwZQZ4QVAA2G6cwATgdhBUCDYTozgNNBWAHQoFqE2vWPa5nODKDuCCsAGtyQDtEax3RmAHVEWAFgir9e1Eldmc4MoA4IKwBMYfez6bnrmc4M4NQIKwBM075l9enMG7PyTa4IgDcirAAwVfXpzKk6UsZ0ZgDVEVYAmOr46cw7cgo1Y9Ems0sC4GUIKwBMd/x05rdXZjCdGUA1hBUAXoHpzABOhLACwGswnRlAbQgrALzGsenMgf5MZwbwP4QVAF6lfcswPTiC6cwA/oewAsDr3Nivjf7Y9X/TmYvLKswuCYCJCCsAvE7ldOaeanl0OvP/+3C9DIPxK0BTRVgB4JWahwTohRvPkp/VooW/Zem1nxi/AjRVhBUAXuvsds310NHl+Gd9sUUrduaZXBEAMxBWAHi1Pw9oq6v6tJLTZWjif37VvvwjZpcEoIERVgB4NYvFohlX9lDXuHDlFZXp9rd/VWkF9w8CmhLCCgCvFxRg08ujUxQR5K/fMg9r+qcbzS4JQAMirADwCQnNg/Xc9X1ksUjvrMrUO6syzC4JQAMhrADwGUM7Ruu+YZ0kSdM+2ajUzMPmFgSgQRBWAPiUO89L1vBulQvG3fH2GuUWlppdEgAPI6wA8CkWi0VPXdNLSdEh2pdfoon/+VUVTpfZZQHwIMIKAJ8TFuivV0anKCTAppVpB/Xk4i1mlwTAgzwaVmbMmKGBAwcqODhYzZo1q3Ubi8VS4zFnzhxPlgWgEWjfMkxPX9tLkjR3abo+/S3L5IoAeIpHw0pZWZmuueYa3XHHHSfdbt68edq3b1/V4+abb/ZkWQAaiYu6x+mO85IlSf/vg3Xaku0wuSIAnuDnyZ0/8sgjkqQ33njjpNs1a9ZMsbGxniwFQCN137BO2rA3X0u35+q2+Wv06cTBigjyN7ssAG7kFWNWJk6cqKioKJ199tmaM2eOXK4TD5YrLS2Vw+Go9gDQdNmsFj03qo9aNQvS7rxi3fNuqlwu7tAMNCamh5XHHntM77//vr755huNGjVKkydP1syZM0+4/axZsxQREVH1SEhIaMBqAXijyJAAvTw6RXY/q77bkqPnvttudkkA3KjeYWX69Om1Doo9/rF69eo67+/BBx/UgAED1Lt3b02ePFmPPvqoZs+efcLtp06dqvz8/KpHZmZmfQ8BQCPUvVWEZlzZQ5L07Dfb9e3m/SZXBMBd6j1mZeLEiRo1atRJt2nXrt3p1qP+/fvL4XBo//79iomJqdFut9tlt9tPe/8AGq+rU1pr3Z7DemvFbt39bqoWThysdlEhZpcF4AzVO6xERUUpKirKE7VIktauXavAwMATTnUGgJN5cERXbcxyaM3uQ7pt/hp9NGGgggM8OpcAgId59G9wRkaGDh48qIyMDDmdTqWmpkqS2rdvr9DQUC1cuFDZ2dkaMGCAgoKCtGTJEj3wwAO69dZbOXsC4LQE+Fn14o1n6dLnf9LW/QW6/4N1ev76PrJYLGaXBuA0WQzD8Niw+TFjxujNN9+s8fySJUt03nnnafHixZo6dap27Nghl8ulpKQkjRs3ThMmTJCfX91ylMPhUEREhPLz8xUeHu7uQwDgo37ZdVDXv7JSFS5DD47oonFDkswuCcBx6vP97dGw0hAIKwBO5K0Vu/TwJxtls1o0f+w5GpjsuUvYAOqnPt/fpk9dBgBPGd2/ra46q5WcLkN/+c9aZR0+YnZJAE4DYQVAo2WxWDTzyh7qGheuvKIy3fH2GpWUO80uC0A9EVYANGqB/ja9PDpFzYL99duefD2ycKPZJQGoJ8IKgEYvoXmwnhvVRxaL9M6qTL2zKsPskgDUA2EFQJNwbsdo3TeskyRp2icblZp52NyCANQZYQVAk3Hnecka3i1GZU6X7nh7jXILS80uCUAdEFYANBkWi0VPXdNLydEh2pdfotvmM+AW8AWEFQBNSligv14e3VfhgX5as/uQ7l6QKqfLp5ebAho9wgqAJqd9y1C98ue+CrBZtXhjth77bJN8fH1MoFEjrABokvontdDT1/aSJL2xfJde+ynd5IoAnAhhBUCTNbJXvB64pIsk6fHPN2vhb1kmVwSgNoQVAE3auCGJGjOwnSRp8nu/6ee0PHMLAlADYQVAk2axWPTQpV11UbdYlTldGv/Wam3fX2B2WQCOQ1gB0OTZrBY9O6q3UtpGylFSoTHzftF+R4nZZQE4irACAKq8h9Crf+6rpKgQ7T18RGPm/aKCknKzywIgwgoAVIkMCdCb/3eOokIDtHmfQ3f++1eVO11mlwU0eYQVADhOQvNgvT7mbAUH2LR0e66mfLieNVgAkxFWAOB3erZuphduOEs2q0Uf/rpHz3y9zeySgCaNsAIAtTi/c0vNuKK7JOm573bonVUZJlcENF2EFQA4gVHntNFdF3SQJD348QZ9t2W/yRUBTRNhBQBO4p4LO+jqlNZyugxN+Pdardtz2OySgCaHsAIAJ2GxWDTrqh4a0iFKR8qd+r83flFGXrHZZQFNCmEFAE7B32bVSzelqGtcuHILy3TzvFU6WFRmdllAk0FYAYA6CLX76Y1bzlarZkFKzy3SuDd/UUm50+yygCaBsAIAddQyPFBv/t/ZCg/0068ZhzVpwVo5XazBAngaYQUA6qF9yzDN/XNfBdis+nLjfj26cCOLxgEeRlgBgHrql9RC/7iulyTpzRW7NXdpmskVAY0bYQUATsOlPeP14IgukqSZi7bo09+yTK4IaLwIKwBwmsYOTtQtg9pJku577zetTMsztyCgkSKsAMBpslgsenBEV13cPVZlTpdufWu1tu0vMLssoNEhrADAGbBZLXrmut7q2zZSjpIKjXl9lfY7SswuC2hUCCsAcIYC/W2a++e+SooOUVZ+icbM+0X5xeVmlwU0GoQVAHCDyJAAvXnLOYoKtWvzPodGv/4zgQVwE8IKALhJQvNgvT3uHDUPCdC6PfkEFsBNCCsA4EadY8P1n/H9CCyAGxFWAMDNCCyAexFWAMADCCyA+xBWAMBDCCyAexBWAMCDCCzAmSOsAICHEViAM+OxsLJr1y6NHTtWiYmJCgoKUnJysqZNm6aysrJq22VkZGjkyJEKCQlRVFSU7rrrrhrbAICvI7AAp89jYWXLli1yuVx6+eWXtXHjRj3zzDOaM2eO/va3v1Vt43Q6NWLECBUVFemnn37SggUL9OGHH2ry5MmeKgsATENgAU6PxTAMo6HebPbs2XrppZeUlpYmSfriiy906aWXKjMzU/Hx8ZKkBQsWaMyYMcrJyVF4ePgp9+lwOBQREaH8/Pw6bQ8AZtuS7dANc3/WwaIy9WgVobfH9lNEsL/ZZQENqj7f3w06ZiU/P1/Nmzev+nnFihXq3r17VVCRpOHDh6u0tFRr1qypdR+lpaVyOBzVHgDgS44/w7J+b75ueo0zLMDJNFhY2blzp55//nndfvvtVc9lZ2crJiam2naRkZEKCAhQdnZ2rfuZNWuWIiIiqh4JCQkerRsAPIHAAtRdvcPK9OnTZbFYTvpYvXp1tddkZWXpoosu0jXXXKNx48ZVa7NYLDXewzCMWp+XpKlTpyo/P7/qkZmZWd9DAACvQGAB6savvi+YOHGiRo0addJt2rVrV/XnrKwsnX/++RowYIBeeeWVatvFxsbq559/rvbcoUOHVF5eXuOMyzF2u112u72+ZQOAVzoWWG6Y+3NVYGEMC1CdRwfY7t27V+eff75SUlL09ttvy2azVWs/NsB2z549iouLkyS9++67uvnmmxlgC6BJYdAtmpr6fH97LKxkZWVp6NChatOmjd56661qQSU2NlZS5dTl3r17KyYmRrNnz9bBgwc1ZswYXXHFFXr++efr9D6EFQCNBYEFTYlXzAb66quvtGPHDn333Xdq3bq14uLiqh7H2Gw2ff755woMDNSgQYN07bXX6oorrtBTTz3lqbIAwGsxhgWoXYOus+IJnFkB0NhwhgVNgVecWQEAnB7OsADVEVYAwAsRWID/IawAgJcisACVCCsA4MV+H1iue2WFsg4fMbssoEERVgDAyx0LLFGhdm3JLtCVLy7Txqx8s8sCGgxhBQB8QOfYcH08YaA6tAzVfkeprp2zQku25phdFtAgCCsA4CNaRwbrgzsGamByCxWVOTXuzdV6e+Vus8sCPI6wAgA+JCLIX2/cco6uTmktp8vQgx9v0KxFm+Vy+fSSWcBJEVYAwMcE+Fk1++qeuvePHSVJL/+Ypr+8s1Yl5U6TKwM8g7ACAD7IYrHorgs66JnresnfZtHn6/fphrkrlVdYanZpgNsRVgDAh13Zp7Xe+r9+Cg/0068Zh3XVS8uVdqDQ7LIAtyKsAICPG5DcQv+9c6BaRwZpd16xrnppuValHzS7LMBtCCsA0Ai0bxmmj+4cpF4JzXS4uFw3vfqzPv0ty+yyALcgrABAIxEdZteC8f01vFuMypwu3fXOWr2wZIcMg5lC8G2EFQBoRIICbHrxxhSNHZwoSZr95VZN/e96lTtdJlcGnD7CCgA0MjarRQ9d2lWPXNZNVou04JdM/d8bv6ighJsgwjcRVgCgkbp5YDu9MrqvgvxtWro9V9fM4SaI8E2EFQBoxC7sGqP3bhug6DBuggjfRVgBgEauR+sIfXTnQHWMOe4miFu4CSJ8B2EFAJqAYzdBHNS+8iaIY9/8hZsgwmcQVgCgiQgP9Ne8MefompTWchniJojwGYQVAGhCAvys+vvVPTX5uJsgTnznV26CCK9GWAGAJsZisegvx90EcdH6bN0wd6VyHCVmlwbUirACAE3UlX1aa/7YfooI8tevGYd1yXNL9dP2XLPLAmogrABAE9Y/qYU+unOgOseGKbewTKNf/1n/+GqrnIxjgRchrABAE5cUHaqPJwzS9ee0kWFIz323Qze+ymUheA/CCgBAgf42zbqqh/45qrdCAmxamXaQy0LwGoQVAECVy3u30qd/GcxlIXgVwgoAoJpkLgvByxBWAAA1nOiy0NLtB8wuDU0QYQUAcEK/vyz059dXcVkIDY6wAgA4KS4LwWyEFQDAKXFZCGYirAAA6ozLQjADYQUAUC9cFkJDI6wAAOqttstCF/9zqT5ft0+GwVkWuBdhBQBw2o6/LJRXVKYJ//lVt7+9hrMscCvCCgDgjCRHh+qTiYN01x/ay89q0Zcb9+vCf/yg91ZncpYFbkFYAQCcMbufTfcO66SFfxmsHq0i5Cip0P0frNPo11Yp82Cx2eXBxxFWAABu0yUuXB/dOVBTL+4su59VP+3I1bBnftS8ZenMGMJp81hY2bVrl8aOHavExEQFBQUpOTlZ06ZNU1lZWbXtLBZLjcecOXM8VRYAwMP8bFbdNjRZi+8+V+ckNteRcqceWbhJ18xZrh05BWaXBx/k56kdb9myRS6XSy+//LLat2+vDRs2aPz48SoqKtJTTz1Vbdt58+bpoosuqvo5IiLCU2UBABpIYlSIFozvr/+sytATX2zRrxmHdck/f9JdF7TXbUOT5W/j5D7qxmI04Oin2bNn66WXXlJaWtr/CrBY9NFHH+mKK644rX06HA5FREQoPz9f4eHhbqoUAOBOWYeP6IGP1mvJ1soVbzvHhmn21b3UozX/OG2q6vP93aCxNj8/X82bN6/x/MSJExUVFaWzzz5bc+bMkcvlasiyAAAeFt8sSK+POVv/HNVbkcH+2pJdoCteXKYnvtiiknKn2eXBy3nsMtDv7dy5U88//7yefvrpas8/9thjuuCCCxQUFKRvv/1WkydPVm5urh588MFa91NaWqrS0tKqnx0Oh0frBgC4h8Vi0eW9W2lw+yhNX7hJC3/L0pwfdurLjdl64qoe6pfUwuwS4aXqfRlo+vTpeuSRR066zS+//KK+fftW/ZyVlaWhQ4dq6NChevXVV0/62qefflqPPvqo8vPz6/X+XAYCAN/y9ab9evDj9drvqPwH6Oj+bXX/RZ0UFuhvcmVoCPW5DFTvsJKbm6vc3NyTbtOuXTsFBgZKqgwq559/vvr166c33nhDVuvJrzwtW7ZMgwcPVnZ2tmJiYmq013ZmJSEhgbACAD7IUVKuWYs2651VmZKk2PBATR7WUVed1Vo2q8Xk6uBJ9Qkr9b4MFBUVpaioqDptu3fvXp1//vlKSUnRvHnzThlUJGnt2rUKDAxUs2bNam232+2y2+31KRkA4KXCA/0166qeGtkzXlP+u14ZB4v11w/W6fVluzT14s46t2O02SXCC3hsNtCxSz9t2rTRW2+9JZvNVtUWGxsrSVq4cKGys7M1YMAABQUFacmSJZo8ebLGjBmjf/7zn3V6H2YDAUDjUFLu1JvLd+lfS3aooKRCkjSkQ5SmXtxFXeP5fG9sPHoZqK7eeOMN3XLLLbW2HXvLxYsXa+rUqdqxY4dcLpeSkpI0btw4TZgwQX5+dTvpQ1gBgMblUFGZnv9uh+av3KVypyGLRbqqT2vdN7yj4iKCzC4PbuIVYaWhEFYAoHHKyCvW37/cos/W7ZMk2f2sGjs4Ubefl6xwBuH6PMIKAKDRWJtxSLMWbdGqXQclSc1DAjTpgg66oV8bVsH1YYQVAECjYhiGvt60X08s3qK0A0WSKpfzv394J13UPVYWCzOHfA1hBQDQKJU7XXr3l0w9+8025RZW3hg3pW2k/nZJZ6W0rblCOrwXYQUA0KgVllbolR92au7SdB05ulz/xd1jdf9FnZUYFWJydagLwgoAoEnY7yjRM19v03urM+UyJD+rRdednaBbz01S2xaEFm9GWAEANClbswv0xBebq+7qbLFIw7vGavy5iVwe8lKEFQBAk7QyLU9zftip74+GFknq06aZxg9J0vBusSzh70UIKwCAJm37/gK9ujRdH63dqzKnS5KU0DxIYwcl6pq+CQqx1/tuM3AzwgoAAJJyCko0f8VuzV+5W4eLyyVJ4YF+urF/W40Z2E4x4YEmV9h0EVYAADjOkTKnPvh1j17/KV3puZXrtPjbLBrZK17jhySpSxzfHw2NsAIAQC2cLkPfbN6vV5em6Zddh6qeH9IhSuOGJOncDlEsMNdACCsAAJxCauZhzV2api/W75Pr6Ddhp5gwjR2SqMt7x8vuZzO3wEaOsAIAQB1lHizWvGW79O4vGSoqq1xgrlmwvy7uHqfLesWrX2JzWZlF5HaEFQAA6in/SLneWZWhN5btUrajpOr52PBAjewVp8t6tVL3VuFcJnITwgoAAKfJ6TK0Mi1Pn6Tu1RcbslVQUlHVlhQVopG94nVZ73glR4eaWKXvI6wAAOAGpRVOfb/1gD79LUvfbNqv0gpXVVv3VuG6vFcrXdorTnERQSZW6ZsIKwAAuFlhaYW+3pStT1KztHR7rpxHR+VaLNI57Zrr8t6tdHH3WEWGBJhcqW8grAAA4EF5haVatCFbn6burTYF2s9q0dCO0bqsd7z+0LmlwgL9TazSuxFWAABoIHsOFeuzdfv0SWqWNu9zVD1vs1rUvVWE+ic1V/+kFurbNpLwchzCCgAAJti+v0Cf/palz9btq1op9xjCS3WEFQAATLb38BH9nJanlWl5Wpl2UBkHi6u1N/XwQlgBAMDLEF6qI6wAAODlsg4f0c/peVq586BWpudpd17N8NK2RbCSokKUFB2qxKgQJUaFKCk6RNGhdp9fnI6wAgCAjzlVeDlemN1PidEhxwWYUCUd/XOI3a8Bqz59hBUAAHzcfkeJduQUKu1AodJyi5R2oEjpuUXac6i46saLtYkJt1cFmMQWIYoMCVCo3abgAD+F2P0UavdTiN2mULufggP8FOBnbbiDOg5hBQCARqq0wqmMvOLjAkxhVZDJKyqr9/4CbFaF2G1VQSY4wHZcqPFTSIBNZyc216U94916HPX5/vaNc0UAAECSZPezqUNMmDrEhNVoyy8uV1puodKPBpldeUVylFSouLRChaUVKiqrUFGpU4WlFSo7euuAMqdLZcUuHSouP+F7lrsMt4eV+iCsAADQSEQE+6tPm0j1aRN5ym3LnS4VlzpVWFahomNhprQyzBQdDTbHnuvVupnniz8JwgoAAE2Qv82qiGCrIoK9f3q0OaNqAAAA6oiwAgAAvBphBQAAeDXCCgAA8GqEFQAA4NUIKwAAwKsRVgAAgFcjrAAAAK9GWAEAAF6NsAIAALwaYQUAAHg1wgoAAPBqhBUAAODVfP6uy4ZhSJIcDofJlQAAgLo69r197Hv8ZHw+rBQUFEiSEhISTK4EAADUV0FBgSIiIk66jcWoS6TxYi6XS1lZWQoLC5PFYnHrvh0OhxISEpSZmanw8HC37hv/Qz83DPq5YdDPDYN+bjie6mvDMFRQUKD4+HhZrScfleLzZ1asVqtat27t0fcIDw/nL0MDoJ8bBv3cMOjnhkE/NxxP9PWpzqgcwwBbAADg1QgrAADAqxFWTsJut2vatGmy2+1ml9Ko0c8Ng35uGPRzw6CfG4439LXPD7AFAACNG2dWAACAVyOsAAAAr0ZYAQAAXo2wAgAAvFqTCisvvviiEhMTFRgYqJSUFC1duvSk2//www9KSUlRYGCgkpKSNGfOnBrbfPjhh+ratavsdru6du2qjz76yFPl+wx39/PcuXM1ZMgQRUZGKjIyUhdeeKFWrVrlyUPwCZ74fT5mwYIFslgsuuKKK9xctW/yRF8fPnxYEyZMUFxcnAIDA9WlSxctWrTIU4fgEzzRz88++6w6deqkoKAgJSQk6J577lFJSYmnDsEn1Kef9+3bpxtuuEGdOnWS1WrV3XffXet2Hv8uNJqIBQsWGP7+/sbcuXONTZs2GZMmTTJCQkKM3bt317p9WlqaERwcbEyaNMnYtGmTMXfuXMPf39/44IMPqrZZvny5YbPZjJkzZxqbN282Zs6cafj5+RkrV65sqMPyOp7o5xtuuMF44YUXjLVr1xqbN282brnlFiMiIsLYs2dPQx2W1/FEPx+za9cuo1WrVsaQIUOMyy+/3MNH4v080delpaVG3759jUsuucT46aefjF27dhlLly41UlNTG+qwvI4n+vntt9827Ha78e9//9tIT083vvzySyMuLs64++67G+qwvE59+zk9Pd246667jDfffNPo3bu3MWnSpBrbNMR3YZMJK+ecc45x++23V3uuc+fOxpQpU2rd/v777zc6d+5c7bnbbrvN6N+/f9XP1157rXHRRRdV22b48OHGqFGj3FS17/FEP/9eRUWFERYWZrz55ptnXrCP8lQ/V1RUGIMGDTJeffVV4+abbyasGJ7p65deeslISkoyysrK3F+wj/JEP0+YMMH4wx/+UG2be++91xg8eLCbqvY99e3n4w0dOrTWsNIQ34VN4jJQWVmZ1qxZo2HDhlV7ftiwYVq+fHmtr1mxYkWN7YcPH67Vq1ervLz8pNucaJ+Nnaf6+feKi4tVXl6u5s2bu6dwH+PJfn700UcVHR2tsWPHur9wH+Spvv700081YMAATZgwQTExMerevbtmzpwpp9PpmQPxcp7q58GDB2vNmjVVl43T0tK0aNEijRgxwgNH4f1Op5/roiG+C33+RoZ1kZubK6fTqZiYmGrPx8TEKDs7u9bXZGdn17p9RUWFcnNzFRcXd8JtTrTPxs5T/fx7U6ZMUatWrXThhRe6r3gf4ql+XrZsmV577TWlpqZ6qnSf46m+TktL03fffacbb7xRixYt0vbt2zVhwgRVVFTo4Ycf9tjxeCtP9fOoUaN04MABDR48WIZhqKKiQnfccYemTJnisWPxZqfTz3XREN+FTSKsHGOxWKr9bBhGjedOtf3vn6/vPpsCT/TzMX//+9/1zjvv6Pvvv1dgYKAbqvVd7uzngoIC3XTTTZo7d66ioqLcX6yPc/fvtMvlUsuWLfXKK6/IZrMpJSVFWVlZmj17dpMMK8e4u5+///57zZgxQy+++KL69eunHTt2aNKkSYqLi9NDDz3k5up9hye+tzz9XdgkwkpUVJRsNluNlJeTk1MjDR4TGxtb6/Z+fn5q0aLFSbc50T4bO0/18zFPPfWUZs6cqW+++UY9e/Z0b/E+xBP9vHHjRu3atUsjR46sane5XJIkPz8/bd26VcnJyW4+Eu/nqd/puLg4+fv7y2azVW3TpUsXZWdnq6ysTAEBAW4+Eu/mqX5+6KGHNHr0aI0bN06S1KNHDxUVFenWW2/VAw88IKu1SYyEqHI6/VwXDfFd2CT+TwUEBCglJUVff/11tee//vprDRw4sNbXDBgwoMb2X331lfr27St/f/+TbnOifTZ2nupnSZo9e7Yee+wxLV68WH379nV/8T7EE/3cuXNnrV+/XqmpqVWPyy67TOeff75SU1OVkJDgsePxZp76nR40aJB27NhRFQgladu2bYqLi2tyQUXyXD8XFxfXCCQ2m01G5eQSNx6Bbzidfq6LBvkudNtQXS93bLrWa6+9ZmzatMm4++67jZCQEGPXrl2GYRjGlClTjNGjR1dtf2xa3D333GNs2rTJeO2112pMi1u2bJlhs9mMJ554wti8ebPxxBNPMHXZA/385JNPGgEBAcYHH3xg7Nu3r+pRUFDQ4MfnLTzRz7/HbKBKnujrjIwMIzQ01Jg4caKxdetW47PPPjNatmxpPP744w1+fN7CE/08bdo0IywszHjnnXeMtLQ046uvvjKSk5ONa6+9tsGPz1vUt58NwzDWrl1rrF271khJSTFuuOEGY+3atcbGjRur2hviu7DJhBXDMIwXXnjBaNu2rREQEGCcddZZxg8//FDVdvPNNxtDhw6ttv33339v9OnTxwgICDDatWtnvPTSSzX2+f777xudOnUy/P39jc6dOxsffvihpw/D67m7n9u2bWtIqvGYNm1aAxyN9/LE7/PxCCv/44m+Xr58udGvXz/DbrcbSUlJxowZM4yKigpPH4pXc3c/l5eXG9OnTzeSk5ONwMBAIyEhwbjzzjuNQ4cONcDReK/69nNtn79t27atto2nvwstRwsBAADwSk1izAoAAPBdhBUAAODVCCsAAMCrEVYAAIBXI6wAAACvRlgBAABejbACAAC8GmEFAAB4NcIKAADwaoQVAADg1QgrAADAqxFWAACAV/v/+aR6PesEwxAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAGdCAYAAADjWSL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0/ElEQVR4nO3de3xU9Z3/8XfI5AIxOeZiEkYDgqWABW+xDUG6QIlcNMStj27WjWZ57LKAi4LhooWytq67BkUFtlIVWX4bStD00bpQ23UjsK0o5R5NK8hCKZGLEILbMCEQc5vz+wNyYCYzgYkzyRnyej4e85A553NOvudIZt58z/d7ToRpmqYAAADCTK/ubgAAAEBnEGIAAEBYIsQAAICwRIgBAABhiRADAADCEiEGAACEJUIMAAAIS4QYAAAQlhzd3YBQcbvdOnHihOLj4xUREdHdzQEAAFfBNE2dPXtWTqdTvXp13NdyzYaYEydOKCMjo7ubAQAAOuHYsWO66aabOqy5ZkNMfHy8pAsnISEhoZtbAwAArkZdXZ0yMjKs7/GOXLMhpu0SUkJCAiEGAIAwczVDQRjYCwAAwhIhBgAAhCVCDAAACEuEGAAAEJYIMQAAICwRYgAAQFgixAAAgLBEiAEAAGGJEAMAAMISIQYAAIQlQgwAAAhLhBgAABCWrtkHQAJd7X+r67T1j1/o/8416boYh+7ql6jM/omKdvBvBQAIBUIM8BVVu77UD9Z/ot/8b027dYl9ovTAHTfqe5k36RvOhKt6KisA4OpEmKZpdncjQqGurk6GYcjlcikhIaG7m4Nr1P9W1+nhVTv1f+eaFNkrQqO/foP6J/fRF/VN2v6nL/RFfZNVOyQ9Xt/LvEkP3HGjboiP6cZWA4B9BfL9TYgBOunzMw16YMVWfVHfpCHp8VpRcJe+lnqdtb6l1a2th77QLyqOa+Onp9TU4pYk9YqQbkrso4E3xKlfUh/1jo5UrCNSsVGRio3qpdioSDl6te+x8dWL46tfx7vMV+dPhI8t6SQCEKgb4mM08paUoO4zkO9vLicBndDqNjXnZ5X6or5JQ/smqGzaCBl9ojxqHJG9NGZwqsYMTpXrfLN+9YcTevuj4/r46Bkd/fN5Hf3z+W5qPQAEx198/Yagh5hAEGKATvj5nmPaVfVnxUVH6vVH7moXYLwZfaL0yIj+emREf50+26g/na7X4dPndOJMgxqaW/Vlc6saW9z6srlVXza71ep2e2zvq7vUVx+q96Kr7Wi9NvtjAYTaN5zde6WDEAMEqKGpVcs2H5Qkzbn36+qfHBfQ9jfEx+iG+BiNGJgciuYBQI/B3E8gQD+vOKZTdY268freKszu393NAYAeixADBOjNnUclSdO+PUAxjshubg0A9FyEGCAAh2rq9b/VZxUVGaHv3nVTdzcHAHo0QgwQgPK9JyVJ93wtRUbvjgfzAgBCixADBOC/91ZLku4b1rebWwIAIMQAV+nM+SbtO1EnSfrO0NRubg0AgBADXKWPjtZKkgbeEKeU63hsAAB0N0IMcJX2fHYhxGT2S+zmlgAAJEIMcNUqjlwIMXffTIgBADsgxABXodVt6vfHz0iS7qInBgBsgRADXIUTZxr0ZbNb0Y5eGnjDdVfeAAAQcgGHmA8++ECTJ0+W0+lURESENmzY4LG+vr5ejz/+uG666Sb17t1bQ4cO1WuvveZR09jYqFmzZiklJUVxcXHKy8vT8ePHPWpqa2tVWFgowzBkGIYKCwt15syZgA8QCIZDp+slSQOS4xTZK6KbWwMAkDoRYs6dO6fbb79dK1as8Ll+zpw5Ki8vV2lpqfbv3685c+Zo1qxZ+uUvf2nVFBUVaf369SorK9PWrVtVX1+v3Nxctba2WjUFBQWqrKxUeXm5ysvLVVlZqcLCwk4cIvDVHT59TtKFmUkAAHsI+CnWkyZN0qRJk/yu3759u6ZMmaIxY8ZIkqZPn66VK1dqz549euCBB+RyubR69WqtXbtWOTk5kqTS0lJlZGRo8+bNmjBhgvbv36/y8nLt2LFDWVlZkqRVq1YpOztbBw4c0ODBgztxqEDn/eliT8wtXEoCANsI+piYUaNG6Z133tHnn38u0zT129/+VgcPHtSECRMkSRUVFWpubtb48eOtbZxOp4YNG6Zt27ZJuhCEDMOwAowkjRgxQoZhWDXeGhsbVVdX5/ECguXwxRBDTwwA2EfQQ8yPf/xj3XrrrbrpppsUHR2tiRMn6tVXX9WoUaMkSdXV1YqOjlZioucMj7S0NFVXV1s1qant74iamppq1XhbvHixNX7GMAxlZGQE+cjQk/3p4uUkemIAwD5CEmJ27Nihd955RxUVFXr55Zc1c+ZMbd68ucPtTNNURMSlAZOX/9lfzeUWLlwol8tlvY4dO/bVDgS46Fxji06fbZQkDaAnBgBsI+AxMR1paGjQD37wA61fv17333+/JOm2225TZWWlXnrpJeXk5Cg9PV1NTU2qra316I2pqanRyJEjJUnp6ek6depUu/2fPn1aaWlpPn92TEyMYmK4FTyC76TrS0lSfKxDCbE8uRoA7CKoPTHNzc1qbm5Wr16eu42MjJTb7ZYkZWZmKioqSps2bbLWnzx5Unv37rVCTHZ2tlwul3bt2mXV7Ny5Uy6Xy6oBukr1xRDT14jt5pYAAC4XcE9MfX29Dh06ZL2vqqpSZWWlkpKS1K9fP40ePVpPPvmkevfurf79+2vLli366U9/qqVLl0qSDMPQ1KlTNW/ePCUnJyspKUnz58/X8OHDrdlKQ4cO1cSJEzVt2jStXLlS0oVZTrm5ucxMQpc76WqQJKUbvbu5JQCAywUcYvbs2aOxY8da7+fOnStJmjJlikpKSlRWVqaFCxfq4Ycf1p///Gf1799fzz33nB599FFrm2XLlsnhcCg/P18NDQ0aN26cSkpKFBkZadWsW7dOs2fPtmYx5eXl+b03DRBKVk9MAj0xAGAnEaZpmt3diFCoq6uTYRhyuVxKSEjo7uYgjP3Thk9UuuOoZn/na5o7np5AAAilQL6/eXYScAW155olSUlx0d3cEgDA5QgxwBX8+VyTJCmREAMAtkKIAa6g9vyFEENPDADYCyEGuIK2EJPYhxADAHZCiAE6YJqmNSaGy0kAYC+EGKAD55pa1dR64UaNSfTEAICtEGKADtReHNQbG9VLvaMjr1ANAOhKhBigA9agXnphAMB2CDFAB9qmV19PiAEA2yHEAB04c75tUC9PrwYAuyHEAB04++WFEBMfQ4gBALshxAAdqG9slSTFxQT8rFQAQIgRYoAOnGtskSRdF8PMJACwG0IM0IH6iyGGnhgAsB9CDNCBc4QYALAtQgzQgXNNbZeTCDEAYDeEGKADDOwFAPsixAAdYGAvANgXIQboAGNiAMC+CDFAB5idBAD2RYgBOnDpchIhBgDshhADdOAcA3sBwLYIMYAfza1uNbW6JUlx0QzsBQC7IcQAfjQ0t1p/jo0ixACA3RBiAD8am93Wn2Mc/KoAgN3wyQz48eXFnpgYRy9FRER0c2sAAN4IMYAfjS0XemK4lAQA9kSIAfy4vCcGAGA/fDoDftATAwD2RogB/Gi82BMTG8WvCQDYEZ/OgB9tPTExDnpiAMCOAg4xH3zwgSZPniyn06mIiAht2LChXc3+/fuVl5cnwzAUHx+vESNG6OjRo9b6xsZGzZo1SykpKYqLi1NeXp6OHz/usY/a2loVFhbKMAwZhqHCwkKdOXMm4AMEOutLemIAwNYC/nQ+d+6cbr/9dq1YscLn+j/96U8aNWqUhgwZovfff1+///3v9fTTTys2NtaqKSoq0vr161VWVqatW7eqvr5eubm5am29dHOxgoICVVZWqry8XOXl5aqsrFRhYWEnDhHoHHpiAMDeAn4gzKRJkzRp0iS/6xctWqT77rtPS5YssZYNHDjQ+rPL5dLq1au1du1a5eTkSJJKS0uVkZGhzZs3a8KECdq/f7/Ky8u1Y8cOZWVlSZJWrVql7OxsHThwQIMHDw602UDA6IkBAHsL6qez2+3Wf/3Xf+nrX/+6JkyYoNTUVGVlZXlccqqoqFBzc7PGjx9vLXM6nRo2bJi2bdsmSdq+fbsMw7ACjCSNGDFChmFYNd4aGxtVV1fn8QK+iraemGimWAOALQX107mmpkb19fV6/vnnNXHiRG3cuFHf/e539eCDD2rLli2SpOrqakVHRysxMdFj27S0NFVXV1s1qamp7fafmppq1XhbvHixNX7GMAxlZGQE89DQAzVffPhjVCQhBgDsKOg9MZL0wAMPaM6cObrjjju0YMEC5ebm6vXXX+9wW9M0PW7t7us27941l1u4cKFcLpf1Onbs2Fc4EkBqcZuSJEcvQgwA2FFQP51TUlLkcDh06623eiwfOnSoNTspPT1dTU1Nqq2t9aipqalRWlqaVXPq1Kl2+z99+rRV4y0mJkYJCQkeL+CraLF6YnhuEgDYUVBDTHR0tL75zW/qwIEDHssPHjyo/v37S5IyMzMVFRWlTZs2WetPnjypvXv3auTIkZKk7OxsuVwu7dq1y6rZuXOnXC6XVQOEWnPrxZ4YQgwA2FLAs5Pq6+t16NAh631VVZUqKyuVlJSkfv366cknn9Rf//Vf6y/+4i80duxYlZeX61e/+pXef/99SZJhGJo6darmzZun5ORkJSUlaf78+Ro+fLg1W2no0KGaOHGipk2bppUrV0qSpk+frtzcXGYmocu0XLw8yuUkALCngEPMnj17NHbsWOv93LlzJUlTpkxRSUmJvvvd7+r111/X4sWLNXv2bA0ePFhvv/22Ro0aZW2zbNkyORwO5efnq6GhQePGjVNJSYkiIy/dj2PdunWaPXu2NYspLy/P771pgFBoudgTw+UkALCnCNM0ze5uRCjU1dXJMAy5XC7Gx6BTnv3Vp/p/v6vSP465Rd+fOKS7mwMAPUIg39/0kwN+tF1OiupFTwwA2BEhBvDj0sBefk0AwI74dAb8aJtiHUlPDADYEiEG8KPtZncM7AUAeyLEAH60PXaAKdYAYE98OgN+MMUaAOyNEAP4Yd3sjoG9AGBLfDoDflizkxjYCwC2RIgB/LDuE0NPDADYEp/OgB9tY2KYYg0A9kSIAfxoeyAHIQYA7IkQA/jhvphiyDAAYE+EGMCPthATEUGKAQA7IsQAfly8Ya96EWIAwJYIMYAfJpeTAMDWCDGAH209MXTEAIA9EWIAP0wxJgYA7IwQA/hx8V53jIkBAJsixAB+MMUaAOyNEAP4YTI7CQBsjRAD+GHdJ6ab2wEA8I0QA/jBze4AwN4IMYAfF68mMSYGAGyKEAP4YY2JIcUAgC0RYgA/mJ0EAPZGiAH8YEwMANgbIQbwo+1md0QYALAnQgxwBdwnBgDsiRAD+HFpTAwhBgDsiBAD+HFpTEw3NwQA4BMhBvDDzWMHAMDWAg4xH3zwgSZPniyn06mIiAht2LDBb+2MGTMUERGh5cuXeyxvbGzUrFmzlJKSori4OOXl5en48eMeNbW1tSosLJRhGDIMQ4WFhTpz5kygzQU6zaQnBgBsLeAQc+7cOd1+++1asWJFh3UbNmzQzp075XQ6260rKirS+vXrVVZWpq1bt6q+vl65ublqbW21agoKClRZWany8nKVl5ersrJShYWFgTYX6DR6YgDA3hyBbjBp0iRNmjSpw5rPP/9cjz/+uN577z3df//9HutcLpdWr16ttWvXKicnR5JUWlqqjIwMbd68WRMmTND+/ftVXl6uHTt2KCsrS5K0atUqZWdn68CBAxo8eHCgzQYCZnKzOwCwtaCPiXG73SosLNSTTz6pb3zjG+3WV1RUqLm5WePHj7eWOZ1ODRs2TNu2bZMkbd++XYZhWAFGkkaMGCHDMKwab42Njaqrq/N4AV9FW08MN7sDAHsKeoh54YUX5HA4NHv2bJ/rq6urFR0drcTERI/laWlpqq6utmpSU1PbbZuammrVeFu8eLE1fsYwDGVkZHzFI0FPx2MHAMDeghpiKioq9G//9m8qKSkJ+F+vpml6bONre++ayy1cuFAul8t6HTt2LLDGA15MxsQAgK0FNcR8+OGHqqmpUb9+/eRwOORwOHTkyBHNmzdPN998syQpPT1dTU1Nqq2t9di2pqZGaWlpVs2pU6fa7f/06dNWjbeYmBglJCR4vICvgvvEAIC9BTXEFBYW6g9/+IMqKyutl9Pp1JNPPqn33ntPkpSZmamoqCht2rTJ2u7kyZPau3evRo4cKUnKzs6Wy+XSrl27rJqdO3fK5XJZNUCo0RMDAPYW8Oyk+vp6HTp0yHpfVVWlyspKJSUlqV+/fkpOTvaoj4qKUnp6ujWjyDAMTZ06VfPmzVNycrKSkpI0f/58DR8+3JqtNHToUE2cOFHTpk3TypUrJUnTp09Xbm4uM5PQZeiJAQB7CzjE7NmzR2PHjrXez507V5I0ZcoUlZSUXNU+li1bJofDofz8fDU0NGjcuHEqKSlRZGSkVbNu3TrNnj3bmsWUl5d3xXvTAMFETwwA2FuE2XYzjGtMXV2dDMOQy+VifAw65Ws/eFctblM7Fo5TuhHb3c0BgB4hkO9vnp0E+MEUawCwN0IM4Efbze5EiAEAWyLEAD5cfpWVMTEAYE+EGMCHy0eKEWIAwJ4IMYAPbo+emG5sCADAL0IM4IP7sp4YHgAJAPZEiAF8uLwnhgwDAPZEiAGugDExAGBPhBjAB8bEAID9EWIAH9zMTgIA2yPEAD4wJgYA7I8QA/hgui/9mZ4YALAnQgzgg0dPTDe2AwDgHyEG8OHyR7vTEwMA9kSIAXxgTAwA2B8hBvChLcRERHDHXgCwK0IM4ENbRwyXkgDAvggxgA9WT0w3twMA4B8hBvDBTU8MANgeIQbwwbxsTAwAwJ4IMYAPjIkBAPsjxAA+tI2J4eGPAGBfhBjAB8bEAID9EWIAH6yb3ZFhAMC2CDGAD4yJAQD7I8QAPpiMiQEA2yPEAD4wJgYA7I8QA/hw6dlJhBgAsCtCDOADU6wBwP4IMYAP1uQkQgwA2BYhBvCB2UkAYH8Bh5gPPvhAkydPltPpVEREhDZs2GCta25u1ve//30NHz5ccXFxcjqd+tu//VudOHHCYx+NjY2aNWuWUlJSFBcXp7y8PB0/ftyjpra2VoWFhTIMQ4ZhqLCwUGfOnOnUQQKBunQ5iRADAHYVcIg5d+6cbr/9dq1YsaLduvPnz+ujjz7S008/rY8++kj/+Z//qYMHDyovL8+jrqioSOvXr1dZWZm2bt2q+vp65ebmqrW11aopKChQZWWlysvLVV5ersrKShUWFnbiEIHAuXkAJADYniPQDSZNmqRJkyb5XGcYhjZt2uSx7JVXXtG3vvUtHT16VP369ZPL5dLq1au1du1a5eTkSJJKS0uVkZGhzZs3a8KECdq/f7/Ky8u1Y8cOZWVlSZJWrVql7OxsHThwQIMHDw602UBAmGINAPYX8jExLpdLERERuv766yVJFRUVam5u1vjx460ap9OpYcOGadu2bZKk7du3yzAMK8BI0ogRI2QYhlXjrbGxUXV1dR4voLNMemIAwPZCGmK+/PJLLViwQAUFBUpISJAkVVdXKzo6WomJiR61aWlpqq6utmpSU1Pb7S81NdWq8bZ48WJr/IxhGMrIyAjy0aAnudgRQ08MANhYyEJMc3OzHnroIbndbr366qtXrDdN0+PGYr5uMuZdc7mFCxfK5XJZr2PHjnW+8ejx3G56YgDA7kISYpqbm5Wfn6+qqipt2rTJ6oWRpPT0dDU1Nam2ttZjm5qaGqWlpVk1p06darff06dPWzXeYmJilJCQ4PECOosxMQBgf0EPMW0B5o9//KM2b96s5ORkj/WZmZmKioryGAB88uRJ7d27VyNHjpQkZWdny+VyadeuXVbNzp075XK5rBoglHgAJADYX8Czk+rr63Xo0CHrfVVVlSorK5WUlCSn06nvfe97+uijj/TrX/9ara2t1hiWpKQkRUdHyzAMTZ06VfPmzVNycrKSkpI0f/58DR8+3JqtNHToUE2cOFHTpk3TypUrJUnTp09Xbm4uM5PQJeiJAQD7CzjE7NmzR2PHjrXez507V5I0ZcoUPfPMM3rnnXckSXfccYfHdr/97W81ZswYSdKyZcvkcDiUn5+vhoYGjRs3TiUlJYqMjLTq161bp9mzZ1uzmPLy8nzemwYIhbb7xAAA7CvCNK/NT+u6ujoZhiGXy8X4GARsy8HTmvL/dunWvgl694lvd3dzAKDHCOT7m2cnAT5Yjx3gNwQAbIuPaMAHk2cnAYDtEWIAH9zuC//1d18iAED3I8QAPlgPgOzmdgAA/CPEAD5cmmLdve0AAPhHiAF8YkwMANgdIQbwgZvdAYD9EWIAH6wxMWQYALAtQgzgAz0xAGB/hBjAB5OeGACwPUIM4INJTwwA2B4hBvCBMTEAYH+EGMAHxsQAgP0RYgAfrAdAkmEAwLYIMYAPPAASAOyPEAP40HY5iQwDAPZFiAF8MK0QQ4oBALsixAA+MCYGAOyPEAP4wJgYALA/QgzgA1OsAcD+CDGAD25rUEz3tgMA4B8hBvCBxw4AgP0RYgAfGNgLAPZHiAF8oCcGAOyPEAP4wAMgAcD+CDGAD8xOAgD7I8QAPlg9Md3cDgCAf4QYoAP0xACAfRFiAB/cF68n9eI3BABsi49owAc3D4AEANsjxAA+cJ8YALC/gEPMBx98oMmTJ8vpdCoiIkIbNmzwWG+app555hk5nU717t1bY8aM0b59+zxqGhsbNWvWLKWkpCguLk55eXk6fvy4R01tba0KCwtlGIYMw1BhYaHOnDkT8AECnWFaA3tJMQBgVwGHmHPnzun222/XihUrfK5fsmSJli5dqhUrVmj37t1KT0/Xvffeq7Nnz1o1RUVFWr9+vcrKyrR161bV19crNzdXra2tVk1BQYEqKytVXl6u8vJyVVZWqrCwsBOHCATu4tUkemIAwMYcgW4wadIkTZo0yec60zS1fPlyLVq0SA8++KAkac2aNUpLS9Obb76pGTNmyOVyafXq1Vq7dq1ycnIkSaWlpcrIyNDmzZs1YcIE7d+/X+Xl5dqxY4eysrIkSatWrVJ2drYOHDigwYMHd/Z4gaty6WZ3pBgAsKugjompqqpSdXW1xo8fby2LiYnR6NGjtW3bNklSRUWFmpubPWqcTqeGDRtm1Wzfvl2GYVgBRpJGjBghwzCsGm+NjY2qq6vzeAGdxc3uAMD+ghpiqqurJUlpaWkey9PS0qx11dXVio6OVmJiYoc1qamp7fafmppq1XhbvHixNX7GMAxlZGR85eNBz8XAXgCwv5DMTvLugjdN84rd8t41vuo72s/ChQvlcrms17FjxzrRcuAC6wGQpBgAsK2ghpj09HRJatdbUlNTY/XOpKenq6mpSbW1tR3WnDp1qt3+T58+3a6Xp01MTIwSEhI8XkBntd3sjggDAPYV1BAzYMAApaena9OmTdaypqYmbdmyRSNHjpQkZWZmKioqyqPm5MmT2rt3r1WTnZ0tl8ulXbt2WTU7d+6Uy+WyaoBQapudxMBeALCvgGcn1dfX69ChQ9b7qqoqVVZWKikpSf369VNRUZGKi4s1aNAgDRo0SMXFxerTp48KCgokSYZhaOrUqZo3b56Sk5OVlJSk+fPna/jw4dZspaFDh2rixImaNm2aVq5cKUmaPn26cnNzmZmELsGYGACwv4BDzJ49ezR27Fjr/dy5cyVJU6ZMUUlJiZ566ik1NDRo5syZqq2tVVZWljZu3Kj4+Hhrm2XLlsnhcCg/P18NDQ0aN26cSkpKFBkZadWsW7dOs2fPtmYx5eXl+b03DRBsJrOTAMD2Isy2W5NeY+rq6mQYhlwuF+NjELAf/nKvfrr9iGZ/52uaO57ePwDoKoF8f/PsJMCHtstJoicGAGyLEAP4cOlmd93bDgCAf4QYwAfGxACA/RFiAB9MZicBgO0RYgAfeAAkANgfIQbwgQdAAoD9EWIAHy71xHRzQwAAfhFiAF+YnQQAtkeIAXy49NgBUgwA2BUhBvDBbd3rjhADAHZFiAF84AGQAGB/hBjAB252BwD2R4gBfGB2EgDYHyEG8MFkTAwA2B4hBvCBMTEAYH+EGMAH7tgLAPZHiAF84AGQAGB/hBjAB2tgr0gxAGBXhBjAh4tXk5idBAA2RogBfGBMDADYHyEG8MEaE8NvCADYFh/RgA88ABIA7I8QA/jgdl/4Lze7AwD7IsQAPlyanQQAsCtCDOBD2+wkLicBgH0RYgAfuNkdANgfIQbwwc0DIAHA9ggxgA88ABIA7I8QA/hATwwA2B8hBvCBMTEAYH+EGMAHk8cOAIDtBT3EtLS06J/+6Z80YMAA9e7dWwMHDtSzzz4rd9vdw3ThX7nPPPOMnE6nevfurTFjxmjfvn0e+2lsbNSsWbOUkpKiuLg45eXl6fjx48FuLuCTdZ8YMgwA2FbQQ8wLL7yg119/XStWrND+/fu1ZMkSvfjii3rllVesmiVLlmjp0qVasWKFdu/erfT0dN177706e/asVVNUVKT169errKxMW7duVX19vXJzc9Xa2hrsJgPt8ABIALA/R7B3uH37dj3wwAO6//77JUk333yz3nrrLe3Zs0fShV6Y5cuXa9GiRXrwwQclSWvWrFFaWprefPNNzZgxQy6XS6tXr9batWuVk5MjSSotLVVGRoY2b96sCRMmBLvZgAeTZycBgO0FvSdm1KhR+p//+R8dPHhQkvT73/9eW7du1X333SdJqqqqUnV1tcaPH29tExMTo9GjR2vbtm2SpIqKCjU3N3vUOJ1ODRs2zKrx1tjYqLq6Oo8X0FlcTgIA+wt6T8z3v/99uVwuDRkyRJGRkWptbdVzzz2nv/mbv5EkVVdXS5LS0tI8tktLS9ORI0esmujoaCUmJraradve2+LFi/XP//zPwT4c9FCmNcW6e9sBAPAv6D0xP/vZz1RaWqo333xTH330kdasWaOXXnpJa9as8ajzvv+GaZpXvCdHRzULFy6Uy+WyXseOHftqB4Ie7dIDIEkxAGBXQe+JefLJJ7VgwQI99NBDkqThw4fryJEjWrx4saZMmaL09HRJF3pb+vbta21XU1Nj9c6kp6erqalJtbW1Hr0xNTU1GjlypM+fGxMTo5iYmGAfDnqotp6YSG4UAwC2FfSemPPnz6tXL8/dRkZGWlOsBwwYoPT0dG3atMla39TUpC1btlgBJTMzU1FRUR41J0+e1N69e/2GGCCYeOwAANhf0HtiJk+erOeee079+vXTN77xDX388cdaunSp/v7v/17ShctIRUVFKi4u1qBBgzRo0CAVFxerT58+KigokCQZhqGpU6dq3rx5Sk5OVlJSkubPn6/hw4dbs5WAUOKxAwBgf0EPMa+88oqefvppzZw5UzU1NXI6nZoxY4Z++MMfWjVPPfWUGhoaNHPmTNXW1iorK0sbN25UfHy8VbNs2TI5HA7l5+eroaFB48aNU0lJiSIjI4PdZKAdemIAwP4izLYbYlxj6urqZBiGXC6XEhISurs5CDP3PP8bfX6mQb987B7dnnF9dzcHAHqMQL6/eXYS4IObm90BgO0RYgAfuNkdANgfIQbwgWcnAYD9EWIAH6xnJ/EbAgC2xUc04AM9MQBgf4QYwAemWAOA/RFiAB/c7raBvaQYALArQgzgg8nlJACwPUIM4AOXkwDA/ggxgA8M7AUA+yPEAD5wszsAsD9CDOADY2IAwP4IMYAPPDsJAOyPEAP4wMBeALA/QgzgQ9vAXu4TAwD2RYgBvLQ9N0miJwYA7IwQA3hxX8owjIkBABsjxABe3B49MYQYALArQgzg5fIQE8FvCADYFh/RgBeTy0kAEBYIMYAXNwN7ASAsEGIALwzsBYDwQIgBvHiMiSHDAIBtEWIAL6b70p/piQEA+yLEAF6YYg0A4YEQA3hhYC8AhAdCDODl8oG9PDsJAOyLEAN4MXmCNQCEBUIM4KWtJ4bxMABgb4QYwEvbmBgyDADYGyEG8HIpxJBiAMDOQhJiPv/8cz3yyCNKTk5Wnz59dMcdd6iiosJab5qmnnnmGTmdTvXu3VtjxozRvn37PPbR2NioWbNmKSUlRXFxccrLy9Px48dD0VzAg2ldTuredgAAOhb0EFNbW6t77rlHUVFR+u///m99+umnevnll3X99ddbNUuWLNHSpUu1YsUK7d69W+np6br33nt19uxZq6aoqEjr169XWVmZtm7dqvr6euXm5qq1tTXYTQY8mIyJAYCw4Aj2Dl944QVlZGToP/7jP6xlN998s/Vn0zS1fPlyLVq0SA8++KAkac2aNUpLS9Obb76pGTNmyOVyafXq1Vq7dq1ycnIkSaWlpcrIyNDmzZs1YcKEYDcbsLit2UmEGACws6D3xLzzzju6++679Vd/9VdKTU3VnXfeqVWrVlnrq6qqVF1drfHjx1vLYmJiNHr0aG3btk2SVFFRoebmZo8ap9OpYcOGWTXeGhsbVVdX5/ECOoOBvQAQHoIeYg4fPqzXXntNgwYN0nvvvadHH31Us2fP1k9/+lNJUnV1tSQpLS3NY7u0tDRrXXV1taKjo5WYmOi3xtvixYtlGIb1ysjICPahoYdgijUAhIeghxi326277rpLxcXFuvPOOzVjxgxNmzZNr732mked98wP0zSvOBuko5qFCxfK5XJZr2PHjn21A0GPxc3uACA8BD3E9O3bV7feeqvHsqFDh+ro0aOSpPT0dElq16NSU1Nj9c6kp6erqalJtbW1fmu8xcTEKCEhweMFdAY9MQAQHoIeYu655x4dOHDAY9nBgwfVv39/SdKAAQOUnp6uTZs2Weubmpq0ZcsWjRw5UpKUmZmpqKgoj5qTJ09q7969Vg0QKtwnBgDCQ9BnJ82ZM0cjR45UcXGx8vPztWvXLr3xxht64403JF34YigqKlJxcbEGDRqkQYMGqbi4WH369FFBQYEkyTAMTZ06VfPmzVNycrKSkpI0f/58DR8+3JqtBISKm8tJABAWgh5ivvnNb2r9+vVauHChnn32WQ0YMEDLly/Xww8/bNU89dRTamho0MyZM1VbW6usrCxt3LhR8fHxVs2yZcvkcDiUn5+vhoYGjRs3TiUlJYqMjAx2kwEP3CcGAMJDhNk2ivEaU1dXJ8Mw5HK5GB+DgPzh+BnlrfidnEasti0c193NAYAeJZDvb56dBHhpG9jLmBgAsDdCDODFGhPDbwcA2Bof04AXk8cOAEBYIMQAXrhPDACEB0IM4MXt5tlJABAOCDGAF3piACA8EGIALzw7CQDCAyEG8EJPDACEB0IM4KWVZycBQFggxABe2gb2RvLbAQC2xsc04KXVCjH8egCAnfEpDXhpaQsxXE0CAFsjxABe2h47EMn0JACwNUIM4OXS5SRCDADYGSEG8EJPDACEB0IM4KWllQdAAkA4IMQAXlrpiQGAsECIAby03SfGQYgBAFsjxABe2qZYczkJAOyNEAN4YWAvAIQHQgzghSnWABAeCDGAF0IMAIQHQgzgxQoxjIkBAFsjxABemGINAOGBEAN4cXM5CQDCAiEG8NLqvvDfXoQYALA1QgzgpdV9IcUwJgYA7I0QA3hhTAwAhAdCDOCl7XISIQYA7I0QA3jhjr0AEB5CHmIWL16siIgIFRUVWctM09Qzzzwjp9Op3r17a8yYMdq3b5/Hdo2NjZo1a5ZSUlIUFxenvLw8HT9+PNTNBdTSSogBgHAQ0hCze/duvfHGG7rttts8li9ZskRLly7VihUrtHv3bqWnp+vee+/V2bNnrZqioiKtX79eZWVl2rp1q+rr65Wbm6vW1tZQNhm41BPDwF4AsLWQhZj6+no9/PDDWrVqlRITE63lpmlq+fLlWrRokR588EENGzZMa9as0fnz5/Xmm29Kklwul1avXq2XX35ZOTk5uvPOO1VaWqpPPvlEmzdvDlWTAUmX7tjLFGsAsLeQhZjHHntM999/v3JycjyWV1VVqbq6WuPHj7eWxcTEaPTo0dq2bZskqaKiQs3NzR41TqdTw4YNs2qAUGmlJwYAwoIjFDstKyvTRx99pN27d7dbV11dLUlKS0vzWJ6WlqYjR45YNdHR0R49OG01bdt7a2xsVGNjo/W+rq7uKx0Deq6Wi9OTHJGEGACws6D3xBw7dkxPPPGESktLFRsb67cuwutfuaZptlvmraOaxYsXyzAM65WRkRF44wFJzRcH9kZHMnkPAOws6J/SFRUVqqmpUWZmphwOhxwOh7Zs2aIf//jHcjgcVg+Md49KTU2NtS49PV1NTU2qra31W+Nt4cKFcrlc1uvYsWPBPjT0EE0Xe2Ki6IkBAFsLeogZN26cPvnkE1VWVlqvu+++Ww8//LAqKys1cOBApaena9OmTdY2TU1N2rJli0aOHClJyszMVFRUlEfNyZMntXfvXqvGW0xMjBISEjxeQGc0t1wMMQ56YgDAzoI+JiY+Pl7Dhg3zWBYXF6fk5GRreVFRkYqLizVo0CANGjRIxcXF6tOnjwoKCiRJhmFo6tSpmjdvnpKTk5WUlKT58+dr+PDh7QYKA8HWbPXEEGIAwM5CMrD3Sp566ik1NDRo5syZqq2tVVZWljZu3Kj4+HirZtmyZXI4HMrPz1dDQ4PGjRunkpISRUZGdkeT0YMwJgYAwkOEaV6cT3qNqaurk2EYcrlcXFpCQPJXbteuqj/rJwV36f7b+nZ3cwCgRwnk+5t/agJemhnYCwBhgRADeLFCDAN7AcDW+JQGvDS3MCYGAMIBn9KAF2YnAUB44FMa8NJ48T4x0VxOAgBb41Ma8HKuqUWSFBfNdH4AsDNCDODlfGOrJCkupltuowQAuEqEGOAyTS1u69lJcdGEGACwM0IMcJnzFy8lSVKfGC4nAYCdEWKAy9Q3Xggx0Y5ezE4CAJvjUxq4zPmmC+NhrmM8DADYHiEGuMz/1TdJkozeUd3cEgDAlRBigMucdDVIkvoasd3cEgDAlRBigMscqqmXJN14fe9ubgkA4Eq48B+gQzX1Kt1xJODtTNP0vbzDbfws72Arf9t09LM62qajFvptX4dt6MT+/O+ucz+rg8X/s/+UJCn7luQOfioAwA4IMQE6caZBJds+6+5mIIScRqxybk3r7mYAAK6AEBOgjKQ+enzs1/yuj4jwv63fVR1s5G9Nxz+ng/35WdXB7jr+WR2t7Mz+/LSkU+e1g+38/Zze0ZG6b3hfJcQysBcA7I4QE6ABKXGaP2FwdzcDAIAej4G9AAAgLBFiAABAWCLEAACAsESIAQAAYYkQAwAAwhIhBgAAhCVCDAAACEuEGAAAEJYIMQAAICwRYgAAQFgixAAAgLBEiAEAAGGJEAMAAMLSNfsUa9M0JUl1dXXd3BIAAHC12r63277HO3LNhpizZ89KkjIyMrq5JQAAIFBnz56VYRgd1kSYVxN1wpDb7daJEycUHx+viIiIoO67rq5OGRkZOnbsmBISEoK6b1zCee4anOeuwXnuGpznrhOqc22aps6ePSun06levToe9XLN9sT06tVLN910U0h/RkJCAr8kXYDz3DU4z12D89w1OM9dJxTn+ko9MG0Y2AsAAMISIQYAAIQlQkwnxMTE6Ec/+pFiYmK6uynXNM5z1+A8dw3Oc9fgPHcdO5zra3ZgLwAAuLbREwMAAMISIQYAAIQlQgwAAAhLhBgAABCWCDGSXn31VQ0YMECxsbHKzMzUhx9+2GH9li1blJmZqdjYWA0cOFCvv/56u5q3335bt956q2JiYnTrrbdq/fr1oWp+2Aj2eV61apW+/e1vKzExUYmJicrJydGuXbtCeQhhIxR/p9uUlZUpIiJCf/mXfxnkVoefUJznM2fO6LHHHlPfvn0VGxuroUOH6t133w3VIYSFUJzn5cuXa/Dgwerdu7cyMjI0Z84cffnll6E6hLAQyHk+efKkCgoKNHjwYPXq1UtFRUU+60L+XWj2cGVlZWZUVJS5atUq89NPPzWfeOIJMy4uzjxy5IjP+sOHD5t9+vQxn3jiCfPTTz81V61aZUZFRZm/+MUvrJpt27aZkZGRZnFxsbl//36zuLjYdDgc5o4dO7rqsGwnFOe5oKDA/MlPfmJ+/PHH5v79+82/+7u/Mw3DMI8fP95Vh2VLoTjXbT777DPzxhtvNL/97W+bDzzwQIiPxN5CcZ4bGxvNu+++27zvvvvMrVu3mp999pn54YcfmpWVlV11WLYTivNcWlpqxsTEmOvWrTOrqqrM9957z+zbt69ZVFTUVYdlO4Ge56qqKnP27NnmmjVrzDvuuMN84okn2tV0xXdhjw8x3/rWt8xHH33UY9mQIUPMBQsW+Kx/6qmnzCFDhngsmzFjhjlixAjrfX5+vjlx4kSPmgkTJpgPPfRQkFodfkJxnr21tLSY8fHx5po1a756g8NYqM51S0uLec8995j//u//bk6ZMqXHh5hQnOfXXnvNHDhwoNnU1BT8BoepUJznxx57zPzOd77jUTN37lxz1KhRQWp1+An0PF9u9OjRPkNMV3wX9ujLSU1NTaqoqND48eM9lo8fP17btm3zuc327dvb1U+YMEF79uxRc3NzhzX+9nmtC9V59nb+/Hk1NzcrKSkpOA0PQ6E8188++6xuuOEGTZ06NfgNDzOhOs/vvPOOsrOz9dhjjyktLU3Dhg1TcXGxWltbQ3MgNheq8zxq1ChVVFRYl58PHz6sd999V/fff38IjsL+OnOer0ZXfBdesw+AvBpffPGFWltblZaW5rE8LS1N1dXVPreprq72Wd/S0qIvvvhCffv29Vvjb5/XulCdZ28LFizQjTfeqJycnOA1PsyE6lz/7ne/0+rVq1VZWRmqpoeVUJ3nw4cP6ze/+Y0efvhhvfvuu/rjH/+oxx57TC0tLfrhD38YsuOxq1Cd54ceekinT5/WqFGjZJqmWlpa9I//+I9asGBByI7Fzjpznq9GV3wX9ugQ0yYiIsLjvWma7ZZdqd57eaD77AlCcZ7bLFmyRG+99Zbef/99xcbGBqG14S2Y5/rs2bN65JFHtGrVKqWkpAS/sWEs2H+n3W63UlNT9cYbbygyMlKZmZk6ceKEXnzxxR4ZYtoE+zy///77eu655/Tqq68qKytLhw4d0hNPPKG+ffvq6aefDnLrw0covrdC/V3Yo0NMSkqKIiMj26XCmpqadumxTXp6us96h8Oh5OTkDmv87fNaF6rz3Oall15ScXGxNm/erNtuuy24jQ8zoTjX+/bt02effabJkydb691utyTJ4XDowIEDuuWWW4J8JPYWqr/Tffv2VVRUlCIjI62aoUOHqrq6Wk1NTYqOjg7ykdhbqM7z008/rcLCQv3DP/yDJGn48OE6d+6cpk+frkWLFqlXr5410qIz5/lqdMV3Yc/6P+UlOjpamZmZ2rRpk8fyTZs2aeTIkT63yc7Oble/ceNG3X333YqKiuqwxt8+r3WhOs+S9OKLL+pf/uVfVF5errvvvjv4jQ8zoTjXQ4YM0SeffKLKykrrlZeXp7Fjx6qyslIZGRkhOx67CtXf6XvuuUeHDh2yQqIkHTx4UH379u1xAUYK3Xk+f/58u6ASGRkp88JklyAeQXjozHm+Gl3yXRi0IcJhqm1a2erVq81PP/3ULCoqMuPi4szPPvvMNE3TXLBggVlYWGjVt03fmzNnjvnpp5+aq1evbjd973e/+50ZGRlpPv/88+b+/fvN559/ninWITjPL7zwghkdHW3+4he/ME+ePGm9zp492+XHZyehONfemJ0UmvN89OhR87rrrjMff/xx88CBA+avf/1rMzU11fzXf/3XLj8+uwjFef7Rj35kxsfHm2+99ZZ5+PBhc+PGjeYtt9xi5ufnd/nx2UWg59k0TfPjjz82P/74YzMzM9MsKCgwP/74Y3Pfvn3W+q74LuzxIcY0TfMnP/mJ2b9/fzM6Otq86667zC1btljrpkyZYo4ePdqj/v333zfvvPNOMzo62rz55pvN1157rd0+f/7zn5uDBw82o6KizCFDhphvv/12qA/D9oJ9nvv3729Kavf60Y9+1AVHY2+h+Dt9OULMBaE4z9u2bTOzsrLMmJgYc+DAgeZzzz1ntrS0hPpQbC3Y57m5udl85plnzFtuucWMjY01MzIyzJkzZ5q1tbVdcDT2Feh59vX5279/f4+aUH8XRlxsCAAAQFjp0WNiAABA+CLEAACAsESIAQAAYYkQAwAAwhIhBgAAhCVCDAAACEuEGAAAEJYIMQAAICwRYgAAQFgixAAAgLBEiAEAAGGJEAMAAMLS/wd95398iWlTxgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameter values\n",
    "p0 = 101325  # pressure\n",
    "T0 = 296  # inlet temperature\n",
    "phi=0.7\n",
    "alpha=0.566\n",
    "mdot_reactants = 5 # kg/m^2/s\n",
    "mdot_products = mdot_reactants  # kg/m^2/s\n",
    "\n",
    "rxnmech = 'gri30.yaml'  # reaction mechanism file\n",
    "gas=ct.Solution(rxnmech)\n",
    "xh2 = alpha / (1 - alpha) * 1  # xh2 = alpha/(1-alpha)*xc\n",
    "fuel = {'CH4': 1, 'H2': xh2}\n",
    "oxidizer = {'O2': 1, 'N2': 3.76}\n",
    "gas.TP = T0, p0\n",
    "gas.set_equivalence_ratio(phi, fuel, oxidizer)  # hold temperature and pressure constant\n",
    "\n",
    "width = 0.1  # m\n",
    "loglevel = 0  # amount of diagnostic output (0 to 5)\n",
    "\n",
    "# Create the flame simulation object\n",
    "flame = ct.CounterflowPremixedFlame(gas=gas, width=width)\n",
    "\n",
    "# Transport model\n",
    "flame.transport_model = 'multicomponent'\n",
    "\n",
    "# Energy equation\n",
    "flame.energy_enabled = True\n",
    "\n",
    "# Set grid refinement parameters\n",
    "flame.set_refine_criteria(ratio=3, slope=0.1, curve=0.2, prune=0.02)\n",
    "\n",
    "# set the boundary flow rates\n",
    "flame.reactants.mdot = mdot_reactants\n",
    "flame.products.mdot = mdot_products\n",
    "\n",
    "flame.set_initial_guess()  # assume adiabatic equilibrium products\n",
    "flame.solve(loglevel, auto=True)\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(flame.grid,flame.velocity)\n",
    "plt.axhline(y=0, color='k', linestyle='--',linewidth=1) #stagnation line\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(flame.grid,flame.T)\n",
    "plt.show()\n",
    "\n",
    "#velocity and temperature\n",
    "#stagnation point(zero); two flows meet\n",
    "#decrease of velocity(use stretch) -> small jump\n",
    "#negative velocity(opposite direction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d58c7b1-5ab8-4905-bc56-3ca524f3246e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
